#=
交换矩阵相关的内容
=#


"""
<rpp| e^{-Δβ*p^2/2m/hbar^2} |rpos> 
= < e^{-Δβm^2(rpp-rpos)^2/(hbar^2Δβ^2*2m*hbar)} > (注: p = mΔx/Δβ)
= exp^{-m(rpp-rpos)^2/2(hbar^2)Δβ}
有时候前面加一个sqrt(m/2πΔβ)
用来归一无相互的Z
λe^2 = 2π*hbar*Δβ/m
exp^{-π(rpp-rpos)^2/λe^2}

free_kinetic(r1, r2)时刻1到时刻2，是增加的
"""
function free_kinetic(rpos::Matrix{T}, rpp::Matrix{T}, vol::Vector{Float64}
    ; λe=_ThermWL) where {T}
    rsqu = map_to_1st((rpp - rpos), vol)
    rsqu = rsqu.^2
    #rsqu1 = @. (rpp - rpos)^2
    return exp(-π*sum(rsqu)/λe/λe)
end


"""
Z = ∫dR ∑_{P} sgn(P) <R| e^{-βH}|P R>
= ∫dR ∑_{P} sgn(P) ∏^{Np} e^{π(r - Pr)^2/λ}
= ∫dR det | e^{π(ri - rj)^2/λ} |
一共Np个粒子,编号(1,2,...,Np)交换他们并且按照逆顺数求和跟求
Np*Np大小的行列式是一样的.其中每个矩阵元ij就是exp^{-π(r_i-r_j)^2/λe^2}
"""
function exch_free_kinetic(rpos::Matrix{T}, rpp::Matrix{T}, vol::Vector{Float64}; λe=_ThermWL) where {T}
    Np, Nd = size(rpos)
    exchmat = zeros(T, Np, Np)
    for rowi in Base.OneTo(Np)
        #rraw = rpp .- rpos[[rowi], :]
        #println("rraw ", rraw)
        #rbnd = map_to_1st.(rraw, vol)
        #println("rbnd ", rbnd)
        #rsqu1 = rbnd.^2
        #rsqu = @. (map_to_1st(rpp - rpos[[rowi], :], vol))^2
        rsqu = map_to_1st((rpp .- rpos[[rowi], :]), vol)
        rsqu = rsqu.^2
        #rsqu1 = @. (rpp - rpos[[rowi], :])^2
        #println(rsqu)
        #println(rsqu1)
        #throw(error())
        rsqu = sum(rsqu, dims=2)
        exchmat[rowi, :] = exp.(-π*rsqu/λe/λe)
    end
    return exchmat
end

